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The application of the diffusion Monte Carlo method to a strongly interacting 
Fermi system as normal liquid 3 He is explored. We show that the fixed-node 
method together with the released-node technique and a systematic method 
to analytically improve the nodal surface constitute an efficient strategy to 
improve the calculation up to a desired accuracy. This methodology shows un- 
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Liquid 3 He has been for many years a benchmark in the field of quantum-many body 
physics. The Fermi statistics of its atoms, combined with the strong correlations induced 
by the hard core of their interatomic potential, has turned it into a paradigm of strongly- 
correlated Fermi systems. At zero temperature, an approximate microscopic description has 
been achieved by means of variational methods, both for the Fermi liquid 3 He as well 
as for its bosonic counterpart liquid 4 He [§J. From a Monte Carlo viewpoint, the quantum 
many-body problem can be tackled in a more ambitious way with the aid of the Green's 
function Monte Carlo (GFMC) and the diffusion Monte Carlo (DMC) methods |§. The 
variational wave function can be used as an input for the Monte Carlo method which, for 
boson systems, is able to solve the Schrodinger equation of the iV-body system providing 
exact results. In Fermi systems as liquid 3 He, the exactness of the method is lost due the 
involved sign problem that makes a straightforward interpretation of the wave function not 
possible. 

The cancellation methods developed up to now to solve this intricate problem have proved 
their efficiency in model problems or with very few particles but become unreliable for real 
many-body systems ||. In the meantime, the approximate fixed-node (FN) J7| method has 
become a standard tool. In the FN-DMC method, the antisymmetry is introduced in the 
trial wave function used for importance sampling imposing its nodal surface as a boundary 
condition. This approach provides upper bounds to the exact eigenvalues, the quality of 
which is related to the accuracy of the nodal surface of the trial wave function. The main 
drawback of the FN-DMC method is the lack of control over the influence of the imposed 
nodal surface on the results obtained, not to say the impossibility of properly correcting 
for such effect. In the present Letter, we come back to this problem using the FN-DMC 
as a main approach but crucially combined with two auxiliary methods: the released-node 
(RN) estimation technique |J and an analytical method able to enhance the quality of any 
given nodal surface. The combination of the above methods provides information on the 
bias due to the imposed nodal surface and a procedure which can evaluate and bring down, 
in principle to any arbitrary required precision, this influence. This complete program has 
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been applied to the study of liquid 3 He bringing, as we will show, the systematic error 
under control and down to levels below the current statistical errors. This has allowed a 
very accurate microscopic calculation of the equation of state of liquid 3 He, including a 
prediction for the negative-pressure region and the spinodal density. 

DMC and GFMC [TI| calculations have provided up to now the best upper bounds 



to the ground-state energy of liquid 3 He . These calculations have unambiguously shown 
the relevance of the Feynman-Cohen-type backflow correlations in the improvement of the 
energy. However, the energy gain appears too small to recover the experimental data and, 
what is more conclusive, the density dependence of the pressure which stresses the curva- 
ture of the equation of state shows clear differences with the experiment. A conclusion that 
naturally emerged from those results was that the nodal surface, originated by backflow cor- 
relations, is not accurate enough and probably new state-dependent correlations ought to be 
considered. Our present results prove that backflow correlations, when properly optimized, 
effectively do provide very accurate nodal surfaces. 

In the DMC method, the Schrodinger equation written in imaginary time is translated 
into a diffusion- like differential equation which can be stochastically solved in an iterative 
procedure. Specific information on the implementation of the DMC method is given in Ref. 
[TT|] . As far as the FN framework [|7| is concerned, the choice of the trial wave function ip(R) 



used for importance sampling is a key point. The simplest model is the Jastrow-Slater wave 
function 

^ = ^ J £> T £> i) (1) 

with if} j = lli<j f{ r ij) a Jastrow wave function and (D\) a Slater determinant of the 
spin-up (spin-down) atoms with single-particle orbitals <p a .(rj) = exp(ik a . • r^). 

In this variational description (ffl), the dynamical correlations induced by the interatomic 
potential are well modelled by the Jastrow factor, and the statistical correlations, implied 
by the antisymmetry, are introduced with a Slater determinant of plane waves which is the 
exact wave function of the free Fermi sea. The two factors account well for the dynamical 



correlations and the Fermi statistics when these effects are independently considered but 
their product is only a relatively poor approximation for a strongly correlated Fermi liquid. 
It is well known, from previous variational and GFMC/DMC calculations PJTOH, that a 
significant improvement on the Jastrow-Slater model is achieved by introducing backflow 
correlations in ip a .(rj), a name which is taken from the Feynman-Cohen famous work on the 
microscopic description of the phonon-roton spectrum in liquid 4 He JL2] . 



At this point, it becomes essential to set up a method for analytically enhancing a given 
model. Such a procedure is already contained in the imaginary-time Schrodinger equation. 
Let us consider a time-dependent wave function 0(R, t), with 0(R, t — 0) = ip(R) the initial 
guess for the trial wave function, satisfying 

-vm-Hw. (2 ) 

A natural choice for a more accurate trial wave function is obtained solving Eq. (Q) at first 
order in t. Near the nodes, which is the relevant region to our purposes, one readily captures 
the main correction to the original i[)a(R) = D^D± in the form 0(R, t) = ^j(R) , 0a(R-) with 
R = R(t) and <9R/ dt = DFj(R), Fj(R) being the drift force coming from the Jastrow 
wave function ipj(R) and D = h 2 /(2m). In this form, the new nodal surface is described 
by the original antisymmetric wave function but with arguments that are shifted due to 
the effect of dynamical correlations. It is worth noting that this approach generates the 
Feynman-Cohen backflow in tp ai (rj) (^ F (R)) as a first order correction to the plane-wave 
orbitals. Recursively, entering with ^ F (R) the next order correction can be analytically 
obtained (see Eq. ||). 

Finally, once a specific model for the nodal surface has been chosen it is necessary to 
establish a method to test its quality. This can be accomplished by means of the released- 
node technique [Q. In the RN approach a superposition of a small boson component in the 
wave function is allowed, with the primary effect of resetting the nodal surface to the exact 
position. This is technically accomplished by introducing a positively-defined guiding wave 
function ^g(R) so th & t the walkers are not confined into a region of definite sign of ^>(R). 



The basic requirements on choosing i/) g (R) are twofold: proximity to |^(R)| away from the 
nodal surface and being positive-defined at the nodes. The choice we have made is 



a being a parameter which controls the crossing frequency. Other different choices can be 
considered but the specific details of i/) g (R) are not relevant since the RN energy is calculated 
projecting out its antisymmetric component. The use of ^g(R) does not introduce any 
systematic bias in the RN energies, which approach the exact eigenvalue when t T — > oo , 
t T being the maximum allowed lifetime after the first crossing. However, the variance of 
the energy grows exponentially with t T due to the boson component, and thus in general 
the asymptotic value cannot be obtained. In contrast, what is straightforwardly available 
is the slope of the energy versus t T at t T — » 0, which provides a direct measure of the 
quality of the input nodal surface (the true antisymmetric ground-state wave function would 
generate a zero slope), and constitutes a means of comparing different trial wave functions. 
In particular, it provides feedback information on whether the next analytical correction to 
?/>a F (R) is necessary. 

We have applied all the above methodology to the study of normal liquid 3 He at zero 
temperature. The results reported have been obtained with N = 66 particles, but we 
have made size checks using also N = 54 and N = 114. In Fermi systems, the kinetic 
energy includes statistical contributions that show an oscillating behavior with N. We have 
observed that this behavior follows very closely that of a discretized Fermi-gas energy, a fact 
that could be expected since such a term appears explicitly in the local kinetic energy. It 
is worth noticing that the case N = 66 is specially well suited for MC calculations as the 
correction amounts only 0.015 K. 

As in previous calculations, we use a short-ranged backflow in the form ip a .(rj) = 
exp(ik ai ■ rf F ), with 



^ g (R) = (^(R) 2 + a 2 ) 1 /2 ? 



(3) 



rf = 



k+3 



jk j 



(4) 
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and r](r) = exp(— ((r — re)/ ^b) 2 )- The two-body correlation factor has been chosen of 
McMillan type, f(r) = exp(-0.5 {b/ rf), and the pairwise HFD-B(HE) Aziz potential |T3 



which has proved high accuracy in liquid 4 He calculations JTTJ, has modelled the atomic 
interactions. At the experimental equilibrium density p[) Xpt = 0.273 cr -3 (er = 2.556 A), we 
have started the calculation with b = 1.15 a and the backflow parameters optimized in Ref. 
jlOl (Ab = 0.14, re = 0.74 a, uj-q = 0.54 a). With this initial set of parameters, the results 
obtained are clearly biased by the trial wave function. Even though the RN approach || 
corrects numerically the shortcomings of ip(R) and, in some applications, allows for an exact 
estimation of the eigenvalue, this is not the case for liquid 3 He. In Fig. 1, the RN energies 
as a function of the released times are shown for the cases Ab = and Ab = 0.14. As one 
can see, at small imaginary times the RN method reveals the presence of corrections to the 
FN energies but a common asymptotic regime is far beyond the scope of the available MC 
data. 

The next step then was looking for the next order correction to backflow correlations, 
as well as for a possibly better set of backflow parameters. We have found that the ones 
we were using correspond to a local minimum of the FN energy, and that a narrower but 
deeper minimum exists with Ab = 0.35 and r-Q and wb unchanged. The resulting energy 
versus released time is also plotted in Fig. 1. The relation of initial slopes, 1 : 0.27 : 0.016 
for Ab = 0, 0.14, 0.35, provides information on the accuracy of ip(R)- In the optimal 
case, Ab = 0.35, the slope is practically inexistent and the energy correction would be 
< 0.01 K if the asymptotic regime could be reached. In order to get additional evidence on 
the size of this correction, and as a closing checkmark of the reliability of our results, we 
have included corrections to the backflow trial wave function using the analytical method 
previously described. It can be shown that these new terms incorporate explicit three-body 
correlations in (p ai (ij) of the form 

rf T = rf + A BT £ V(r jk ) & - , (5) 

with Ti = J2i^iV{ r u) r u- We have carried out a FN-DMC calculation with this new trial 
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wave function at P(j xpt and the result for the energy correction has been found < 0.01 K. 
Both this analytical check and the numerical findings provided by the RN method point out 
the excellent description that backflow correlations make of the nodal surface in liquid 3 He. 

The FN energies with Ab = (no backflow), Ab = 0.14, and Ab = 0.35 are reported in 
Table I, together with the corresponding kinetic energy obtained as the difference between 



the total energy and a pure estimation [14| of the potential energy. The comparison with 
the experimental energy , also contained in the table, shows the successive improvement 
of the FN-DMC result until an excellent agreement with Ab = 0.35. Concerning the kinetic 
energy, a sizeable difference between theory and experiment [[16] survives, a fact that has been 
generally attributed to long-range wings in the high-g inelastic response that are difficult to 
incorporate effectively in the experimental analysis. 

The FN-DMC calculation has been extended to a wide range of densities ranging from 
the spinodal point up to a maximum value p = 0.403 cr~ 3 , located near to the experimental 
freezing density pj xpt = 0.394 a~ 3 . We have used AT = 114 only at the highest density 
and below that N = 66 have proved to be accurate enough. Among the three variational 
parameters entering in the backflow wave function (|4|) only Ab shows a density dependence 
which is nearly linear in the range studied (Ab = 0.42 at p = 0.403 o --3 ). The results are 
displayed in Fig. 2 in comparison with the experimental data of Ref. [[HJ. The solid line 
in the same figure is a third-degree polynomial fit to our data with \ 2 l v = 1-2. According 
to this fit, the equilibrium density is po = 0.274(1) o~~ 3 and the energy at this density 
(E/N)o = —2.464(7) K, in close agreement with experimental data. 

The quality of the equation of state is even more stressed by looking at its derivatives. 
In Fig. 3, the behavior of the pressure and the sound velocity with the density is shown 
in comparison with experimental data from Refs. |i^,|l7|, |l8| . The theoretical prediction 
for both quantities, derived form the polynomial fit to E/N(p) (Fig. 1), shows again an 
excellent agreement with the experimental data from equilibrium up to freezing. The sound 
velocity that at Po xpt is c = 182.2(6) m/ sec, in close agreement with the experimental 
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value c expt = 182.9 m/ sec ||17|| , goes down to zero at the spinodal point. The location of 
this point has been previously obtained both from extrapolation of experimental data at 
positive pressures [jnj and from density-functional theories ||20|| . The present microscopic 
calculation allows for an accurate calculation, free from extrapolation uncertainties, that 
locates the spinodal point at a density p s = 0.202(2) <r~ 3 corresponding to a negative pressure 
P s = —3.09(20) atm, much closer to the equilibrium than in liquid 4 He where P s = —9.30(15) 
atm m. 

In conclusion, we have analyzed the possibilities of the diffusion Monte Carlo method 
in the study of Fermi systems. The FN-DMC method, combined with the RN technique 
to evaluate the effect of the nodal surface model used in the trial function, and with a 
systematic method to analytically improve the nodal surface, constitutes a closed loop able 
to improve the quality of the antisymmetric wave function up to a required precision. The 
fact that this general approach has allowed to deal with a strongly interacting Fermi liquid 
suggests that it could be also useful for tackling other Fermi systems. We have applied this 
methodology to the study of liquid 3 He. The effect of corrections beyond the backflow terms 
has been evaluated and found to be less than 0.01 K. The equation of state so obtained 
presents an accuracy comparable to the result obtained in bosonic liquid 4 He. The precision 
proved by the method allows for posterior studies as the characterization of spin-polarized 
liquid 3 He. Preliminary calculations for the fully-polarized phase at p[j xpt indicate a less 
binded system with an energy E/N = —2.22(4) K. 

This research has been partially supported by DGES (Spain) Grant N° PB96-0170-C03- 
02. We also acknowledge the supercomputer facilities provided by the CEPBA. 
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TABLES 

TABLE I. FN-DMC total and kinetic energies for liquid 3 He at function of the 

backflow parameter A B - The experimental values for the total and kinetic energies are taken form 
Ref. [15 1 and Ref. [jn|, respectively. 



E/N (K) T/N (K) 



A B = -2.128 ±0.015 12.603 ±0.031 

A B = 0.14 -2.330 ±0.014 12.395 ± 0.035 

A B = 0.35 -2.477 ±0.014 12.239 ± 0.030 

Expt -2.473 8.1 ± 1.7 
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FIGURES 

FIG. 1. Relesead-node energies as a function of the released time. Open circles, triangles and 
full circles correspond to Ab = (no backflow), Ab = 0.14, and Ab = 0.35, respectively. The lines 
are linear fits to the MC data. 




12 



FIG. 2. Energy per particle of normal liquid 3 He as a function of the density. The full circles 
are the FN-DMC results (the error bars are depicted only when larger than the size of the symbol), 



and the open circles are experimental data from Ref. [15|. The line is a polynomial fit to the MC 
data. 
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FIG. 3. Pressure and sound velocity as a function of the density. The lines are the FN-DMC 
results and the circles, triangles and squares are experimental data from Refs. IB], [17], and [18], 
respectively. 
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